// Define_Extensions
*estimates benchmark specifications: Tables 2-4

cd "$path_name"

// PRELIMINARIES
	# delimit ;
	set more off;
	clear all;
	matrix drop _all;

	set mem 800m;

	cap log close _all;

	*identify dates for first-stage estimation;
	cap ma drop fsdates;
	global fsdates "year>=1994 & year<=1996";

	*load data;
	global data_file="CE_data.dta";
	use "$data_file";

	* Use as switch for showing elasticities
	global display_beta 0;

	*variable names for goods;
	global variables shelt fdhome vehpch oppub health fdaway util equpmt adultcl shoesoth furnsh
	tvrdot feeard housop  chldrn educa persca alcbev cashco tobacc;

	*variable names for goods with suffix r indicating deflated expenditure;
	global variablesr sheltr fdhomer vehpchr oppubr utilr healthr equpmtr fdawayr  tvrdotr adultclr 
	feeardr cashcor furnshr educar shoesothr housopr  alcbevr chldrnr  perscar tobaccr ;

	*variable names for nondurable goods;
	global variables_nd shelt_nd fdhome_nd vehpch_nd oppub_nd health_nd fdaway_nd util_nd equpmt_nd adultcl_nd shoesoth_nd furnsh_nd
	tvrdot_nd feeard_nd housop_nd  chldrn_nd educa_nd persca_nd alcbev_nd cashco_nd tobacc_nd;

	*variable names for nondurable goods with suffix r;
	global variablesr_nd sheltr_nd fdhomer_nd vehpchr_nd oppubr_nd healthr_nd fdawayr_nd utilr_nd equpmtr_nd adultclr_nd shoesothr_nd furnshr_nd
	tvrdotr_nd feeardr_nd housopr_nd  chldrnr_nd educar_nd perscar_nd alcbevr_nd cashcor_nd tobaccr_nd;


	*year and income class dummies;
	gen d95=year==1995;
	gen dyr96=year==1996;
	gen dyr06=year==2006;
	gen dyr07=year==2007;
	xi i.inclass, pre(D);

	*calculate nominal expenditure share of each good in CE for the first-stage years 1994-1996;
	*Reported in table 1;
	foreach var in $variables totexp {;
		qui gen nom_m_`var'=.;
		foreach num of numlist 1994/1996 {;
			quietly sum `var' [aw=nweight] if year==`num';
			quietly replace nom_m_`var'=r(mean) if year==`num';
		};
	};
	egen tag_year=tag(year);
	foreach var in $variables {;
		qui gen CEshare_`var'r= nom_m_`var'/nom_m_totexp if tag_year==1;
		qui egen temp=mean(CEshare_`var'r);
		qui replace CEshare_`var'r=temp;
		drop temp;
	};
	drop tag_year;


	sort id;

// FIRST STAGE
 
	**********************************************************************************;
	*First stage: estimate expenditure elasticities using 1994-1996 ($fsdates) data;
	*IV expenditure with income category and log income;
	capture program drop first_stage;

	program define first_stage, eclass;

	cap matrix drop Beta;
	cap matrix drop SE;
	cap matrix drop FStage;

	foreach var in $variablesr {;
		quietly ivreg dev_`var' (ltotexpr = Dinclass* lfincatar) age2d age3d npers2d npers3d nearn2d dyr* [pw=nweight]   if $fsdates; 
		matrix b=e(b);
		matrix Beta=[nullmat(Beta) \ b[1,"ltotexpr".."nearn2d"]];
		if 1 == 1 {;
			matrix se=e(V);
			matrix se=se["ltotexpr","ltotexpr"];
			matrix se=sqrt(se[1,1]);
			matrix SE=[nullmat(SE) \ se];
			qui sum CEshare_`var';
			matrix CEshare=[nullmat(CEshare) \ r(mean)*100];
		};
	};
	matrix rownames Beta= $variablesr;
	cap matrix FStage=[CEshare, Beta[1...,"ltotexpr"], SE];
	cap matrix rownames FStage= $variablesr;
	cap matrix colnames FStage= Share Beta SE;
	ereturn matrix Beta=Beta;
	cap ereturn matrix SE=SE;
	cap ereturn matrix FStage=FStage;

	end;

	**********************************************************************************;

	**********************************************************************************;
	*First stage with alternative IV: estimate expenditure elasticities using 1994-1996 ($fsdates) data;
	*IV Q45 expenditure with Q23 expenditure;

	capture program drop first_stage_qtr;

	program define first_stage_qtr, eclass;

	cap matrix drop Beta;
	cap matrix drop SE;

	foreach var in $variablesr {;
		quietly ivreg dev_Q45_`var' (ltotexpr_Q45 = ltotexpr_Q23) age2d age3d npers2d npers3d nearn2d dyr* [pw=nweight]   if $fsdates; 
		matrix b=e(b);
		matrix Beta=[nullmat(Beta) \ b[1,"ltotexpr_Q45".."nearn2d"]];
		matrix se=e(V);
		matrix se=se["ltotexpr_Q45","ltotexpr_Q45"];
		matrix se=sqrt(se[1,1]);
		matrix SE=[nullmat(SE) \ se];
	};


	matrix FStage_qtr=[Beta[1...,"ltotexpr_Q45"], SE];

	matrix rownames Beta= $variablesr;
	matrix rownames SE= $variablesr;
	matrix rownames FStage_qtr= $variablesr;
	matrix colnames FStage_qtr= BetaQ SEQ;

	ereturn matrix Beta=Beta;
	ereturn matrix SE=SE;
	ereturn matrix FStage_qtr=FStage_qtr;

	end;

	**********************************************************************************;


	**********************************************************************************;
	*First stage with only nondurable expenditures;
	*IV expenditure with income category and log income;
	capture program drop first_stage_nd;

	program define first_stage_nd, eclass;

	cap matrix drop Beta;
	cap matrix drop SE;
	cap matrix drop FStage;

	foreach var in $variablesr {;
		quietly ivreg dev_`var'_nd (ltotexpr_nd = Dinclass* lfincatar) age2d age3d npers2d npers3d nearn2d dyr* [pw=nweight]   if $fsdates; 
		matrix b=e(b);
		matrix Beta=[nullmat(Beta) \ b[1,"ltotexpr_nd".."nearn2d"]];
		*THIS IS CURRENTLY A PATCH FOR THE FACT THAT $DISPLAY_BETA IS WEIRD.;
		if 1==1 {;
			matrix se=e(V);
			matrix se=se["ltotexpr_nd","ltotexpr_nd"];
			matrix se=sqrt(se[1,1]);
			matrix SE=[nullmat(SE) \ se];
			qui sum CEshare_`var';
			matrix CEshare=[nullmat(CEshare) \ r(mean)*100];
		};
	};
	matrix rownames Beta= $variablesr;
	matrix colnames Beta= ltotexpr age2d age3d npers2d npers3d nearn2d ;
	cap matrix FStage=[CEshare, Beta[1...,"ltotexpr"], SE];
	cap matrix rownames FStage= $variablesr;
	cap matrix colnames FStage= Share Beta SE;
	ereturn matrix Beta=Beta;
	cap ereturn matrix SE=SE;
	cap ereturn matrix FStage_nd=FStage;

	end;

// SECOND STAGE - create dataset

	****************************************************************************
	capture program drop second_stage_dataset;

	program define second_stage_dataset, rclass;
	
	preserve;
		
	//di `"`goodlim'"'
	*replace id for the random bootstrap draw;
	qui replace id=_n;

	*estimate first stage beta's ;
	first_stage;
	matrix BetaAn=e(Beta);

	/* Here I will run ther three possible extensions to (3)
		3: Run whole test only using variables with a "Reasonable" number of zeros
		3': Run log (3) Dropping zeros
		3'': Run log (3) 0 = 0 + 1
		
		Note that when I say I have "log" I really mean log deviation from the mean.
		3' vars are l_var, 3'' are lzero_var
		The criterion for a "reasonable number of zeros" is <3000. */

	// Create global for the arguments.
	global lvariablesr;
	foreach var in $variablesr {;
		cap drop l_`var';
		cap drop mean;
		qui egen mean = wtmean(sheltr), weight(nweight) by(year);
		qui gen l_`var' = log(`var') - log(mean);
		qui gen lzero_`var' = l_`var';
		qui replace lzero_`var' = 1 - log(mean) if lzero_`var'==.;
		global lvariablesr $lvariablesr l_`var';
	};

	pause on;

	*three-year windows CHANGED - I make more windows!;
	qui gen decade=1 if year<=1982 & year>=1980;
	qui replace decade=2 if year<=1985 & year>=1983;
	qui replace decade=3 if year>=1986 & year<=1990;
	qui replace decade=4 if year>=1991 & year<=1993;
	qui replace decade=5 if year>=1995 & year<=1997;
	qui replace decade=6 if year>=1998 & year<=2000;
	qui replace decade=7 if year>=2001 & year<=2004;
	qui replace decade=8 if year>=2005 & year<=2007;
	qui replace decade=9 if year>=2008;
	qui drop if decade==.;

	// Save for summary statistics
	save "SummaryData.dta", replace;
	
	*adjust dependent variable for fitted demographics;
	*generate beta (expenditure elasticity) variable for each good using first-stage results;
	foreach var in $variablesr {;
		cap matrix drop bdemo;
		matrix bdemo=BetaAn["`var'","age2d".."nearn2d"];
		cap drop demo;
		matrix score demo=bdemo;
		qui gen xAn`var'=dev_`var'-demo;
		qui gen l_xAn`var'=l_`var'-demo;
		qui gen lzero_xAn`var'=lzero_`var'-demo;
		cap drop demo;
		matrix betaAn=BetaAn["`var'","ltotexpr"];
		qui gen betaAn`var'=betaAn[1,1];
	};


	keep xAn* l_xAn* lzero_xAn* CEshare_* betaAn* year decade inclass id nweight nsh_mean*;

	* dev_*r;
	*drop dev_totexpr;
	*drop dev_Q*;

	*reshape long xAn xQtr xND betaAn betaQtr betaND nsh_mean nsh_meanND dev_, i(id) j(goods) str;
	drop nsh_meanND*;
	reshape long xAn l_xAn lzero_xAn betaAn nsh_mean, i(id) j(goods) str;
/*	
	// Don't know why it need be a global but hey i'll take it:
	
	global GL `"`goodlim'"';
	keep $GL ;
*/
	*generate good-year dummies;
	sort goods;
	qui egen goodyear=group(goods year);
	
	*generate decadeXincome-class dummies
	char inclass[omit] 1;
	char decade[omit] 1;
	qui xi i.decade*i.inclass, noomit pre(D);

	*DdecX_x_y will be dummy for decade x and income y;
	foreach num_dec of numlist 2/9 {;
		foreach num_inc of numlist 2/5 {;
			rename DdecX~_`num_dec'_`num_inc' DdecX_`num_dec'_`num_inc';
		};
	};

	// b_x_y will be beta interaction for decade x and income y;
	foreach dec of numlist 1/9 {;
		foreach num of numlist 2/5  {;
			qui gen bAn_`dec'_`num'=betaAn*Dinclass_`num'*Ddecade_`dec';
		};
		drop DdecXinc_`dec'_1;
	};
	
	// Save Data File
	save "CleanShape.dta", replace;
	
	restore;
	
	end;
	
// SECOND STAGE - run regressions.
	
	****************************************************************************
	capture program drop second_stage_reg;

	program define second_stage_reg, rclass;

	use "CleanShape.dta", replace;
	
	// Define the limiter, and calculate consumption share
	local i = 1;
	local share = 0;
	local goodlim = "if";
	if "`1'" == "" {;
		di "Well there are no vars here.";
		local share = 1;
	};
	else{;
		while "``i''" != ""{;
			local share = `share' + CEshare_``i'';
			local g = "``i''";
			local goodlim = `"`goodlim' goods == "``i''" |"';
			local i = `i' + 1;
		};
		// Add something false to not mess up the "Or"
		local goodlim = `"`goodlim' 0==1"';
		
		// Drop vars - Not sure why global needed, but okay.
		global GL `"`goodlim'"';
		keep $GL ;
	};

	// Return consumption share of the chosen goods
	return scalar Cshare=`share';
	
	*OLS;
	qui areg xAn DdecX* bAn_* [aw=nweight], absorb(goodyear); 
	matrix OLS=e(b);
	return scalar N_OLS = e(N);

	qui areg l_xAn DdecX* bAn_* [aw=nweight], absorb(goodyear); 
	matrix l_OLS=e(b);
	return scalar N_l_OLS = e(N);

	qui areg lzero_xAn DdecX* bAn_* [aw=nweight], absorb(goodyear); 
	matrix lzero_OLS=e(b);
	return scalar N_lzero_OLS = e(N);

		foreach num of numlist 1/9 {;
			foreach num2 of numlist 3 5 {;
				matrix temp=OLS[1,"bAn_`num'_`num2'"];	
				return scalar OLSb_`num'_`num2'=temp[1,1];
				cap matrix drop temp;
				
				matrix temp=l_OLS[1,"bAn_`num'_`num2'"];	
				return scalar l_OLSb_`num'_`num2'=temp[1,1];
				cap matrix drop temp;
				
				matrix temp=lzero_OLS[1,"bAn_`num'_`num2'"];	
				return scalar lzero_OLSb_`num'_`num2'=temp[1,1];
				cap matrix drop temp;

			};
		};

	end;
	
// Define program which starts by constructing dataset and then regresses;	
	
	capture program drop second_stage_all;

	program define second_stage_all, rclass;
	
	preserve;
	//di `"`goodlim'"'
	*replace id for the random bootstrap draw;
	qui replace id=_n;

	*estimate first stage beta's ;
	first_stage;
	matrix BetaAn=e(Beta);

	/* Here I will run ther three possible extensions to (3)
		3: Run whole test only using variables with a "Reasonable" number of zeros
		3': Run log (3) Dropping zeros
		3'': Run log (3) 0 = 0 + 1
		
		Note that when I say I have "log" I really mean log deviation from the mean.
		3' vars are l_var, 3'' are lzero_var
		The criterion for a "reasonable number of zeros" is <3000. */

	// Create global for the arguments.
	global lvariablesr;
	foreach var in $variablesr {;
		cap drop l_`var';
		cap drop mean;
		qui egen mean = wtmean(sheltr), weight(nweight) by(year);
		qui gen l_`var' = log(`var') - log(mean);
		qui gen lzero_`var' = l_`var';
		qui replace lzero_`var' = 1 - log(mean) if lzero_`var'==.;
		global lvariablesr $lvariablesr l_`var';
	};
	
	pause on;

	*three-year windows CHANGED - I make more windows!;
	qui gen decade=1 if year<=1982 & year>=1980;
	qui replace decade=2 if year<=1985 & year>=1983;
	qui replace decade=3 if year>=1986 & year<=1990;
	qui replace decade=4 if year>=1991 & year<=1993;
	qui replace decade=5 if year>=1995 & year<=1997;
	qui replace decade=6 if year>=1998 & year<=2000;
	qui replace decade=7 if year>=2001 & year<=2004;
	qui replace decade=8 if year>=2005 & year<=2007;
	qui replace decade=9 if year>=2008;
	qui drop if decade==.;

	// Save for summary statistics
	save "SummaryData.dta", replace;
	
	*adjust dependent variable for fitted demographics;
	*generate beta (expenditure elasticity) variable for each good using first-stage results;
	foreach var in $variablesr {;
		cap matrix drop bdemo;
		matrix bdemo=BetaAn["`var'","age2d".."nearn2d"];
		cap drop demo;
		matrix score demo=bdemo;
		qui gen xAn`var'=dev_`var'-demo;
		qui gen l_xAn`var'=l_`var'-demo;
		qui gen lzero_xAn`var'=lzero_`var'-demo;
		cap drop demo;
		matrix betaAn=BetaAn["`var'","ltotexpr"];
		qui gen betaAn`var'=betaAn[1,1];
	};

	// Include (Residualized) Income-based measurements of inequality
	reg fincatar age2d age3d npers2d npers3d nearn2d;
	predict inc_resid, residuals;
	qui su fincatar;
	local avg_inc = r(mean);
	replace inc_resid = `avg_inc' + inc_resid;
	
	foreach dec of numlist 1/9{;
		qui su inc_resid [aw=nweight] if decade==`dec' & inclass==5;
		local avg_resid_5_`dec' = r(mean);
		qui su inc_resid [aw=nweight] if decade==`dec' & inclass==3;
		local avg_resid_3_`dec' = r(mean);
		qui su inc_resid [aw=nweight] if decade==`dec' & inclass==1;
		local avg_resid_1_`dec' = r(mean);
		return scalar resid_Hi_v_Low_Dec`dec' = `avg_resid_5_`dec''/`avg_resid_1_`dec'';
		return scalar resid_Med_v_Low_Dec`dec' = `avg_resid_3_`dec''/`avg_resid_1_`dec'';
		return scalar resid_Hi_v_Med_Dec`dec' = `avg_resid_5_`dec''/`avg_resid_3_`dec'';
		
		
		qui su fincatar [aw=nweight] if decade==`dec' & inclass==5;
		local avg_inc_5_`dec' = r(mean);
		qui su fincatar [aw=nweight] if decade==`dec' & inclass==3;
		local avg_inc_3_`dec' = r(mean);
		qui su fincatar [aw=nweight] if decade==`dec' & inclass==1;
		local avg_inc_1_`dec' = r(mean);
		return scalar inc_Hi_v_Low_Dec`dec' = `avg_inc_5_`dec''/`avg_inc_1_`dec'';
		return scalar inc_Med_v_Low_Dec`dec' = `avg_inc_3_`dec''/`avg_inc_1_`dec'';
		return scalar inc_Hi_v_Med_Dec`dec' = `avg_inc_5_`dec''/`avg_inc_3_`dec'';			
	};
	
	// Report Differences
	// Just 1
		return scalar logch_resid_Hi_v_Low_1 = 
			log(`avg_resid_5_1'/`avg_resid_1_1');
		return scalar logch_resid_Hi_v_Med_1 = 
			log(`avg_resid_5_1'/`avg_resid_3_1');	
		return scalar logch_resid_Med_v_Low_1 = 
			log(`avg_resid_3_1'/`avg_resid_1_1');

	
	// 1 to 4
		return scalar logch_resid_Hi_v_Low_4_to_1 = 
			log(`avg_resid_5_4'/`avg_resid_1_4') - log(`avg_resid_5_1'/`avg_resid_1_1');
		return scalar logch_resid_Hi_v_Med_4_to_1 = 
			log(`avg_resid_5_4'/`avg_resid_3_4') - log(`avg_resid_5_1'/`avg_resid_3_1');	
		return scalar logch_resid_Med_v_Low_4_to_1 = 
			log(`avg_resid_3_4'/`avg_resid_1_4') - log(`avg_resid_3_1'/`avg_resid_1_1');

		return scalar logch_inc_Hi_v_Low_4_to_1 = 
			log(`avg_inc_5_4'/`avg_inc_1_4') - log(`avg_inc_5_1'/`avg_inc_1_1');
		return scalar logch_inc_Hi_v_Med_4_to_1 = 
			log(`avg_inc_5_4'/`avg_inc_3_4') - log(`avg_inc_5_1'/`avg_inc_3_1');	
		return scalar logch_inc_Med_v_Low_4_to_1 = 
			log(`avg_inc_3_4'/`avg_inc_1_4') - log(`avg_inc_3_1'/`avg_inc_1_1');	
			
	// 4 to 6
		return scalar logch_resid_Hi_v_Low_6_to_4 = 
			log(`avg_resid_5_6'/`avg_resid_1_6') - log(`avg_resid_5_4'/`avg_resid_1_4');
		return scalar logch_resid_Hi_v_Med_6_to_4 = 
			log(`avg_resid_5_6'/`avg_resid_3_6') - log(`avg_resid_5_4'/`avg_resid_3_4');	
		return scalar logch_resid_Med_v_Low_6_to_4 = 
			log(`avg_resid_3_6'/`avg_resid_1_6') - log(`avg_resid_3_4'/`avg_resid_1_4');

		return scalar logch_inc_Hi_v_Low_6_to_4 = 
			log(`avg_inc_5_6'/`avg_inc_1_6') - log(`avg_inc_5_4'/`avg_inc_1_4');
		return scalar logch_inc_Hi_v_Med_6_to_4 = 
			log(`avg_inc_5_6'/`avg_inc_3_6') - log(`avg_inc_5_4'/`avg_inc_3_4');	
		return scalar logch_inc_Med_v_Low_6_to_4 = 
			log(`avg_inc_3_6'/`avg_inc_1_6') - log(`avg_inc_3_4'/`avg_inc_1_4');	
			
	// 6 to 8
		return scalar logch_resid_Hi_v_Low_8_to_6 = 
			log(`avg_resid_5_8'/`avg_resid_1_8') - log(`avg_resid_5_6'/`avg_resid_1_6');
		return scalar logch_resid_Hi_v_Med_8_to_6 = 
			log(`avg_resid_5_8'/`avg_resid_3_8') - log(`avg_resid_5_6'/`avg_resid_3_6');	
		return scalar logch_resid_Med_v_Low_8_to_6 = 
			log(`avg_resid_3_8'/`avg_resid_1_8') - log(`avg_resid_3_6'/`avg_resid_1_6');

		return scalar logch_inc_Hi_v_Low_8_to_6 = 
			log(`avg_inc_5_8'/`avg_inc_1_8') - log(`avg_inc_5_6'/`avg_inc_1_6');
		return scalar logch_inc_Hi_v_Med_8_to_6 = 
			log(`avg_inc_5_8'/`avg_inc_3_8') - log(`avg_inc_5_6'/`avg_inc_3_6');	
		return scalar logch_inc_Med_v_Low_8_to_6 = 
			log(`avg_inc_3_8'/`avg_inc_1_8') - log(`avg_inc_3_6'/`avg_inc_1_6');	
			
	// 8 to 9
		return scalar logch_resid_Hi_v_Low_9_to_8 = 
			log(`avg_resid_5_9'/`avg_resid_1_9') - log(`avg_resid_5_8'/`avg_resid_1_8');
		return scalar logch_resid_Hi_v_Med_9_to_8 = 
			log(`avg_resid_5_9'/`avg_resid_3_9') - log(`avg_resid_5_8'/`avg_resid_3_8');	
		return scalar logch_resid_Med_v_Low_9_to_8 = 
			log(`avg_resid_3_9'/`avg_resid_1_9') - log(`avg_resid_3_8'/`avg_resid_1_8');

		return scalar logch_inc_Hi_v_Low_9_to_8 = 
			log(`avg_inc_5_9'/`avg_inc_1_9') - log(`avg_inc_5_8'/`avg_inc_1_8');
		return scalar logch_inc_Hi_v_Med_9_to_8 = 
			log(`avg_inc_5_9'/`avg_inc_3_9') - log(`avg_inc_5_8'/`avg_inc_3_8');	
		return scalar logch_inc_Med_v_Low_9_to_8 = 
			log(`avg_inc_3_9'/`avg_inc_1_9') - log(`avg_inc_3_8'/`avg_inc_1_8');
			
	// 1 to 8
		
		return scalar logch_resid_Hi_v_Low_8_to_1 = 
			log(`avg_resid_5_8'/`avg_resid_1_8') - log(`avg_resid_5_1'/`avg_resid_1_1');
		return scalar logch_resid_Hi_v_Med_8_to_1 = 
			log(`avg_resid_5_8'/`avg_resid_3_8') - log(`avg_resid_5_1'/`avg_resid_3_1');	
		return scalar logch_resid_Med_v_Low_8_to_1 = 
			log(`avg_resid_3_8'/`avg_resid_1_8') - log(`avg_resid_3_1'/`avg_resid_1_1');

		return scalar logch_inc_Hi_v_Low_8_to_1 = 
			log(`avg_inc_5_8'/`avg_inc_1_8') - log(`avg_inc_5_1'/`avg_inc_1_1');
		return scalar logch_inc_Hi_v_Med_8_to_1 = 
			log(`avg_inc_5_8'/`avg_inc_3_8') - log(`avg_inc_5_1'/`avg_inc_3_1');	
		return scalar logch_inc_Med_v_Low_8_to_1 = 
			log(`avg_inc_3_8'/`avg_inc_1_8') - log(`avg_inc_3_1'/`avg_inc_1_1');	
			
	// 1 to 9
		return scalar logch_resid_Hi_v_Low_9_to_1 = 
			log(`avg_resid_5_9'/`avg_resid_1_9') - log(`avg_resid_5_1'/`avg_resid_1_1');
		return scalar logch_resid_Hi_v_Med_9_to_1 = 
			log(`avg_resid_5_9'/`avg_resid_3_9') - log(`avg_resid_5_1'/`avg_resid_3_1');	
		return scalar logch_resid_Med_v_Low_9_to_1 = 
			log(`avg_resid_3_9'/`avg_resid_1_9') - log(`avg_resid_3_1'/`avg_resid_1_1');

		return scalar logch_inc_Hi_v_Low_9_to_1 = 
			log(`avg_inc_5_9'/`avg_inc_1_9') - log(`avg_inc_5_1'/`avg_inc_1_1');
		return scalar logch_inc_Hi_v_Med_9_to_1 = 
			log(`avg_inc_5_9'/`avg_inc_3_9') - log(`avg_inc_5_1'/`avg_inc_3_1');	
		return scalar logch_inc_Med_v_Low_9_to_1 = 
			log(`avg_inc_3_9'/`avg_inc_1_9') - log(`avg_inc_3_1'/`avg_inc_1_1');
	
	keep xAn* l_xAn* lzero_xAn* CEshare_* betaAn* year decade inclass id nweight nsh_mean*;
		
	* dev_*r;
	*drop dev_totexpr;
	*drop dev_Q*;

	*reshape long xAn xQtr xND betaAn betaQtr betaND nsh_mean nsh_meanND dev_, i(id) j(goods) str;
	drop nsh_meanND*;
	reshape long xAn l_xAn lzero_xAn betaAn nsh_mean, i(id) j(goods) str;
/*	
	// Don't know why it need be a global but hey i'll take it:
	
	global GL `"`goodlim'"';
	keep $GL ;
*/
	*generate good-year dummies;
	sort goods;
	qui egen goodyear=group(goods year);
	
	*generate decadeXincome-class dummies
	char inclass[omit] 1;
	char decade[omit] 1;
	qui xi i.decade*i.inclass, pre(D);

	*DdecX_x_y will be dummy for decade x and income y;
	foreach num_inc of numlist 2/5 {;
		qui gen DdecX_1_`num_inc'=Dinclass_`num_inc';
	};
	foreach num_dec of numlist 2/9 {;
		foreach num_inc of numlist 2/5 {;
			rename DdecX~_`num_dec'_`num_inc' DdecX_`num_dec'_`num_inc';
		};
	};
	//b_x_y will be beta interaction for decade x and income y;
	//Decade 1;
	foreach num of numlist 2/5 {;
		qui gen bAn_1_`num'=betaAn*Dinclass_`num';
	};
	//Decades 2-9
	/*foreach d of numlist 2/9 {;
		foreach num of numlist 2/5  {;
			di "Yo: `d'";
			di "Ddecade_`d'";
			qui gen bAn_`d'_`num'=betaAn*Dinclass_`num' * Ddecade_`d';
		}; 
		drop DdecXinc_`dec'_1;
	};*/
	di "indeed";
	foreach d of numlist 2/9 {;
		foreach num of numlist 2/5 {;
			qui gen bAn_`d'_`num'=betaAn*Dinclass_`num' * Ddecade_`d';
		};
	};


	// Save Data File
	save "CleanShape.dta", replace;
	
	*** The Second Part Starts Here ***

	use "CleanShape.dta", replace;
	// Somehow... CleanShape is not getting used???
	// Define the limiter, and calculate consumption share
	local i = 1;
	local share = 0;
	local goodlim = "if";
	if "`1'" == "" {;
		di "Well there are no vars here.";
		local share = 1;
	};
	else{;
		while "``i''" != ""{;
			local share = `share' + CEshare_``i'';
			local g = "``i''";
			local goodlim = `"`goodlim' goods == "``i''" |"';
			local i = `i' + 1;
		};
		// Add something false to not mess up the "Or"
		local goodlim = `"`goodlim' 0==1"';

		// Drop vars - Not sure why global needed, but okay.
		global GL `"`goodlim'"';
		use "CleanShape.dta", replace;
		keep $GL;
	};
	
	// Return consumption share of the chosen goods
	return scalar Cshare=`share';
	
	*OLS;
	qui areg xAn DdecX* bAn_* [aw=nweight], absorb(goodyear); 
	matrix OLS=e(b);
	return scalar N_OLS = e(N);
	
	qui areg l_xAn DdecX* bAn_* [aw=nweight], absorb(goodyear); 
	matrix l_OLS=e(b);
	return scalar N_l_OLS = e(N);

	qui areg lzero_xAn DdecX* bAn_* [aw=nweight], absorb(goodyear); 
	matrix lzero_OLS=e(b);
	return scalar N_lzero_OLS = e(N);

		foreach num of numlist 1/9 {;
			foreach num2 of numlist 3 5 {;
				matrix temp=OLS[1,"bAn_`num'_`num2'"];	
				return scalar OLSb_`num'_`num2'=temp[1,1];
				cap matrix drop temp;
				
				matrix temp=l_OLS[1,"bAn_`num'_`num2'"];	
				return scalar l_OLSb_`num'_`num2'=temp[1,1];
				cap matrix drop temp;
				
				matrix temp=lzero_OLS[1,"bAn_`num'_`num2'"];	
				return scalar lzero_OLSb_`num'_`num2'=temp[1,1];
				cap matrix drop temp;

				matrix temp=OLS[1,"DdecX_`num'_`num2'"];	
				return scalar OLSD_`num'_`num2'=temp[1,1];
				cap matrix drop temp;
			};
		};
		
	// Also Generate some quick differences to include in the formatted tables;
	// (Between 1991-93/98-00, 1998-00/2005-07, 2005-07/2008-10);
		foreach num of numlist 3 5 {;
			matrix temp_i = OLS[1,"bAn_4_`num'"];
			matrix temp_f = OLS[1,"bAn_6_`num'"];
			return scalar IneqChange_`num'_4to6 = temp_f[1,1] - temp_i[1,1] ; 
			matrix drop temp_i temp_f;
			
			matrix temp_i = OLS[1,"bAn_6_`num'"];
			matrix temp_f = OLS[1,"bAn_8_`num'"];
			return scalar IneqChange_`num'_6to8 = temp_f[1,1] - temp_i[1,1] ; 
			matrix drop temp_i temp_f;
			
			matrix temp_i = OLS[1,"bAn_8_`num'"];
			matrix temp_f = OLS[1,"bAn_9_`num'"];
			return scalar IneqChange_`num'_8to9 = temp_f[1,1] - temp_i[1,1] ; 
			matrix drop temp_i temp_f;
			
			matrix temp_i = OLS[1,"bAn_1_`num'"];
			matrix temp_f = OLS[1,"bAn_9_`num'"];
			return scalar IneqChange_`num'_1to8 = temp_f[1,1] - temp_i[1,1] ; 
			matrix drop temp_i temp_f;
		
			matrix temp_i = OLS[1,"bAn_1_`num'"];
			matrix temp_f = OLS[1,"bAn_8_`num'"];
			return scalar IneqChange_`num'_1to8 = temp_f[1,1] - temp_i[1,1] ; 
			matrix drop temp_i temp_f;
			
			matrix temp_i = OLS[1,"DdecX_4_`num'"];
			matrix temp_f = OLS[1,"DdecX_6_`num'"];
			return scalar MismesChange_`num'_4to6 = temp_f[1,1] - temp_i[1,1] ; 
			matrix drop temp_i temp_f;
			
			matrix temp_i = OLS[1,"DdecX_6_`num'"];
			matrix temp_f = OLS[1,"DdecX_8_`num'"];
			return scalar MismesChange_`num'_6to8 = temp_f[1,1] - temp_i[1,1] ; 
			matrix drop temp_i temp_f;
			
			matrix temp_i = OLS[1,"DdecX_8_`num'"];
			matrix temp_f = OLS[1,"DdecX_9_`num'"];
			return scalar MismesChange_`num'_8to9 = temp_f[1,1] - temp_i[1,1] ; 
			matrix drop temp_i temp_f;
		
			matrix temp_i = OLS[1,"DdecX_1_`num'"];
			matrix temp_f = OLS[1,"DdecX_8_`num'"];
			return scalar MismesChange_`num'_1to8 = temp_f[1,1] - temp_i[1,1] ; 
			matrix drop temp_i temp_f;
		
			matrix temp_i = OLS[1,"DdecX_1_`num'"];
			matrix temp_f = OLS[1,"DdecX_9_`num'"];
			return scalar MismesChange_`num'_1to9 = temp_f[1,1] - temp_i[1,1] ; 
			matrix drop temp_i temp_f;
		};
	
	matrix temp_i = OLS[1,"bAn_1_3"];
	matrix temp_f = OLS[1,"bAn_1_5"];
	return scalar PanB1 = temp_f[1,1] - temp_i[1,1]; 
	matrix drop temp_i temp_f;	
	
	matrix temp_i = OLS[1,"bAn_4_3"];
	matrix temp_f = OLS[1,"bAn_4_5"];
	return scalar PanB4 = temp_f[1,1] - temp_i[1,1]; 
	matrix drop temp_i temp_f;	
	
	matrix temp_i1 = OLS[1,"bAn_4_5"];
	matrix temp_i2 = OLS[1,"bAn_4_3"];
	matrix temp_f1 = OLS[1,"bAn_6_5"];
	matrix temp_f2 = OLS[1,"bAn_6_3"];
	return scalar PanB4to6 = (temp_f1[1,1] - temp_i1[1,1]) - (temp_f2[1,1] - temp_i2[1,1]); 
	matrix drop temp_i1 temp_f1 temp_i2 temp_f2;	
	
	matrix temp_i1 = OLS[1,"bAn_6_5"];
	matrix temp_i2 = OLS[1,"bAn_6_3"];
	matrix temp_f1 = OLS[1,"bAn_8_5"];
	matrix temp_f2 = OLS[1,"bAn_8_3"];
	return scalar PanB6to8 = (temp_f1[1,1] - temp_i1[1,1]) - (temp_f2[1,1] - temp_i2[1,1]); 
	matrix drop temp_i1 temp_f1 temp_i2 temp_f2;	
	
	matrix temp_i1 = OLS[1,"bAn_8_5"];
	matrix temp_i2 = OLS[1,"bAn_8_3"];
	matrix temp_f1 = OLS[1,"bAn_9_5"];
	matrix temp_f2 = OLS[1,"bAn_9_3"];
	return scalar PanB8to9 = (temp_f1[1,1] - temp_i1[1,1]) - (temp_f2[1,1] - temp_i2[1,1]); 
	matrix drop temp_i1 temp_f1 temp_i2 temp_f2;	
	
	matrix temp_1 = OLS[1,"bAn_8_5"];
	matrix temp_2 = OLS[1,"bAn_8_3"];
	return scalar PanB1to8 = (temp_1[1,1] - temp_2[1,1]); 
	matrix drop temp_1 temp_2;		
	
	matrix temp_1 = OLS[1,"bAn_9_5"];
	matrix temp_2 = OLS[1,"bAn_9_3"];
	return scalar PanB1to9 = (temp_1[1,1] - temp_2[1,1]); 
	matrix drop temp_1 temp_2;	

	matrix temp_i = OLS[1,"DdecX_1_3"];
	matrix temp_f = OLS[1,"DdecX_1_5"];
	return scalar PanB1_mismes = temp_f[1,1] - temp_i[1,1]; 
	matrix drop temp_i temp_f;	
	
	matrix temp_i = OLS[1,"DdecX_4_3"];
	matrix temp_f = OLS[1,"DdecX_4_5"];
	return scalar PanB4_mismes = temp_f[1,1] - temp_i[1,1]; 
	matrix drop temp_i temp_f;	
	
	matrix temp_i1 = OLS[1,"DdecX_4_5"];
	matrix temp_i2 = OLS[1,"DdecX_4_3"];
	matrix temp_f1 = OLS[1,"DdecX_6_5"];
	matrix temp_f2 = OLS[1,"DdecX_6_3"];
	return scalar PanB4to6_mismes = (temp_f1[1,1] - temp_i1[1,1]) - (temp_f2[1,1] - temp_i2[1,1]); 
	matrix drop temp_i1 temp_f1 temp_i2 temp_f2;	
	
	matrix temp_i1 = OLS[1,"DdecX_6_5"];
	matrix temp_i2 = OLS[1,"DdecX_6_3"];
	matrix temp_f1 = OLS[1,"DdecX_8_5"];
	matrix temp_f2 = OLS[1,"DdecX_8_3"];
	return scalar PanB6to8_mismes = (temp_f1[1,1] - temp_i1[1,1]) - (temp_f2[1,1] - temp_i2[1,1]); 
	matrix drop temp_i1 temp_f1 temp_i2 temp_f2;	
	
	matrix temp_i1 = OLS[1,"DdecX_8_5"];
	matrix temp_i2 = OLS[1,"DdecX_8_3"];
	matrix temp_f1 = OLS[1,"DdecX_9_5"];
	matrix temp_f2 = OLS[1,"DdecX_9_3"];
	return scalar PanB8to9_mismes = (temp_f1[1,1] - temp_i1[1,1]) - (temp_f2[1,1] - temp_i2[1,1]); 
	matrix drop temp_i1 temp_f1 temp_i2 temp_f2;	
	
	matrix temp_1 = OLS[1,"DdecX_8_5"];
	matrix temp_2 = OLS[1,"DdecX_8_3"];
	return scalar PanB1to8_mismes = (temp_1[1,1] - temp_2[1,1]); 
	matrix drop temp_1 temp_2;
	
	matrix temp_1 = OLS[1,"DdecX_9_5"];
	matrix temp_2 = OLS[1,"DdecX_9_3"];
	return scalar PanB1to9_mismes = (temp_1[1,1] - temp_2[1,1]); 
	matrix drop temp_1 temp_2;	

	restore;

	end;
	
